% Based on replication codes for Dario Caldara and Christophe Kamps (2017), The Analytics of SVARs: A
% Unified Framework to Measure Fiscal Multipliers, Review of Economic Studies (2017) 84, 1015–1040
%
% Copyright: 2017 Dario Caldara and Christophe Kamps
% Copyright: 2020-2023 Benjamin Born, Francesco D'Ascanio, Gernot J. Mueller, Johannes Pfeifer

clear all
tic
%-----------------
% Housekeeping
%----------------
if ~isfolder('results')
    mkdir('.','results');
end
if ~isfolder('Figures')
    mkdir('.','Figures');
end


options = optimoptions('fsolve','Display','none');
%--------------------
% Model Specification
%--------------------
p               = 4;      % Number of lags
nex_            = 1;      % Constant
nP              = 2; % Number of fiscal variables
T0              = p;      % Size of pre-sample for Minnesota Prior
Horizon         = 41;     % Horizon for Impulse Responses
nd              = 20000;  % Number of draws in MC chain
bburn           = 0*nd; % Share of draws to burn; 0 here due to conjugate prior not inducing a Markov Chain

ptileVEC = [0.05 0.50 0.95]; % Percentiles for posterior distributions

do_BP        = 1;   % Execute Blanchard Perotti (2002) identification
do_PF        = 1;   % Execute penalty function identification

% Number of periods for penalty function identification
nperPF = 1;


%--------------------
% Load Data
%--------------------

splits = 1; % 0: no split, 1: euro country-time obs only
[data] = load_panel_data_for_VAR_FE_intra_euro_FX(p,splits);
YY=data.VAR.Y;
XXact=data.VAR.X;
setup.var_names=data.var_names;
n               = size(setup.var_names,2); % Number of variables
GDPRATIO =[1 0.2 0.2];


% Fiscal elasticities for Blanchard Perotti identification
BP_EETA_SIMPLE      = zeros(n-2,nP);
BP_EETA_SIMPLE(1,1) = 1.7;
BP_EETA_GENERAL     = [1.7 0; 1.25 -0.50; 0 0; 0 0];

%-------------------------
% Generate Minnesota Prior
%-------------------------

% tau    : Overall tightness
% d      : Scaling down the variance for the coefficients of a distant lag
% w      : Number of observations used for obtaining the prior for the
%          covariance matrix of error terms (usually fixed to 1).
% lambda : Tuning parameter for coefficients for constant.
% mu     : Tuning parameter for the covariance between coefficients*/

%--------------------
% Dummy Observations
%--------------------

tau	   =   0.2;
d	   =   2.5;
w	   =   0;
lambda =   0.5;
mu	   =   0.5;

%get "presample" moments
YY0     =   YY(1:T0,:);
ybar    =   mean(YY0)';
sbar    =   std(YY0)';
pre_sample_moments  =   [ybar sbar];

%------------------------------------------
% Generate matrices with dummy observations
%------------------------------------------

hyp = [tau; d; w; lambda; mu];
[YYdum, XXdum] = varprior_h(n,p,nex_,hyp,pre_sample_moments);

%--------------------
% Build matrix of actual observations
%--------------------
YYact = YY;
nobs    = size(YY,1);  % number of observations
XXact = [XXact ones(nobs,1)];

%-------------------------
% Estimation Preliminaries: get OLS estimates and construct companion form
%-------------------------

J = [eye(n);repmat(zeros(n),p-1,1)]; % Page 12 RWZ; % J: shock linear combination selector matrix
F = zeros(n*p,n*p);    % Matrix for Companion Form
I  = eye(n);
for lag=1:p-1
    F(lag*n+1:(lag+1)*n,(lag-1)*n+1:lag*n) = I;
end

X = [XXdum; XXact];
Y = [YYdum; YYact];
T = size(X, 1);
ndum = size(XXdum, 1);

% Compute OLS estimates
B = (X'*X)\(X'*Y);       % Point estimates
U = Y-X*B;               % Residuals
Sigmau = U'*U/(T-p*n-1); % Covariance matrix of residuals
F(1:n,1:n*p)    = B(1:n*p,:)';


[F21Ord1BP, IRF_scaled_BP, structural_shocks_BP]= vm_mult_simple(F,J,Horizon,BP_EETA_SIMPLE,U,T,n,p,nP,GDPRATIO); 
structural_shocks.BP_T=NaN(data.linear_indices.original_dimensions);
structural_shocks.BP_G=NaN(data.linear_indices.original_dimensions);
structural_shocks.BP_T(data.linear_indices.linear_indices)=structural_shocks_BP(end-nobs+1:end,1);
structural_shocks.BP_G(data.linear_indices.linear_indices)=structural_shocks_BP(end-nobs+1:end,2);

[~, A0PF, Ltilde] = penalty_fun_simple(Sigmau,n,p,F,J,Horizon,nperPF);

PF_EETA_SIMPLE      = zeros(n-2,nP);
pos.T=1;
pos.G=2;
pos.Y=3;
pos.FG=6;
shock.BC=1;
shock.G=2;
shock.T=3;

% get elasticity of T and G w.r.t. other variables (Y and Forecast here)
PF_EETA_SIMPLE(pos.Y-nP,pos.T) = -A0PF(pos.Y,shock.T)/A0PF(pos.T,shock.T); % First row tax, third row GDP; Second column is tax shock;
PF_EETA_SIMPLE(pos.Y-nP,pos.G) = -A0PF(pos.Y,shock.G)/A0PF(pos.G,shock.G); % Second row spending, third row GDP; Third column is spending shock;
PF_EETA_SIMPLE(pos.FG-nP,pos.T) = -A0PF(pos.FG,shock.T)/A0PF(pos.T,shock.T); % Second row spending, third row GDP; Third column is spending shock;
PF_EETA_SIMPLE(pos.FG-nP,pos.G) = -A0PF(pos.FG,shock.G)/A0PF(pos.G,shock.G); % Second row spending, third row GDP; Third column is spending shock;
[F21Ord1PF, IRF_scaled_PF,structural_shocks_PF] = vm_mult_simple(F,J,Horizon,PF_EETA_SIMPLE,U,T,n,p,nP,GDPRATIO);

structural_shocks.PF_T=NaN(data.linear_indices.original_dimensions);
structural_shocks.PF_G=NaN(data.linear_indices.original_dimensions);
structural_shocks.PF_T(data.linear_indices.linear_indices)=structural_shocks_PF(end-nobs+1:end,1);
structural_shocks.PF_G(data.linear_indices.linear_indices)=structural_shocks_PF(end-nobs+1:end,2);

corr(structural_shocks.PF_T(:),structural_shocks.BP_T(:),"rows","pairwise")
corr(structural_shocks.PF_G(:),structural_shocks.BP_G(:),"rows","pairwise")
save('results/structural_shocks_point','structural_shocks');